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ABSTRACT 

We propose a fast and efficient bispectrum statistic for Cosmic Microwave Background 
(CMB) temperature anisotropies to constrain the amplitude of the primordial non- 
Gaussian signal measured in terms of the non-linear coupling parameter /nl- We 
show how the method can achieve a remarkable computational advantage by focussing 
on subsets of the multipole configurations, where the non-Gaussian signal is more 
concentrated. The detection power of the test, increases roughly linearly with the 
maximum multipole, as shown in the ideal case of an experiment without noise and 
gaps. The CPU-time scales as instead of for the full bispectrum which 
for Planck resolution (ma-x ~ 3000 means an improvement in speed of a factor 10^ 
compared to the full bispectrum analysis with minor loss in precision. We find that 
the introduction of a galactic cut partially destroys the optimality of the configuration, 
which will then need to be dealt with in the future. We find for an ideal experiment 
with £niax = 2000 that upper limits of /nl < 8 can be obtained at Icr. For the case 
of the WMAP experiment, we would be able to put limits of |/nl| < 40 if no galactic 
cut were present. Using the real data with galactic cut, we obtain an estimate of 
—80 < /nl < 80 and —160 < /nl < 160 at 1 and 2o- respectively. 

Key words: cosmic microwave background - cosmology: theory - methods: numerical 
- methods: statistical - cosmology: observations 


1 INTRODUCTION 

In recent years a number of papers have focussed on the statistical nature of the fluctuations in the Cosmic Microwave Back¬ 
ground radiation (CMB). In the standard inflationary scenario the quantum fluctuations of the inflaton scalar field follow a 
nearly Gaussian distribution, with small deviations arise by considering second-order terms of the equations JAcquaviva et al. 2003| 
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IMaldacena 2003ll : slow-roll conditions necessarily entail that deviations from Gaussianity is very low. Nonetheless, the sub¬ 
sequent gravitational evolution unavoidably enhances the primordial non-Gaussian signal up to the largest scales where the 
primordial seeds were produced, giving rise to a non-linearity parameter (see below) /nl ~ Cl(l). The situation changes in 
the presence of a second scalar field during inflation; under these circumstances it has been shown by (Bartolo Matarrese and 
Riotto 2002), that non-Gaussianity can be transferred from the isocurvature to the adiabatic mode, leading to non-negligible 
values for /nl (see also IIRernardea.u and TIza.nl [Rigopoulos, Shellard and van Tent 2005| |Seery and Lidsey 2005| l). Alterna¬ 
tive scenarios for the production of the primordial seeds, such as the curvaton{see lIMollerach 1990| [Engvist and Sloth 2002| 
|Lyth and Wands 2002|lMoroi and Takahashi 2001l[Lvth. Ungarelli and Wands 2003^ ) and the inhomogeneous reheating mech¬ 
anisms (see Dvali, Gruzinov and Zaldarriaga 2004) may also lead to higher values of /nl- A general review for the primor¬ 
dial non-Gaussian scenarios can be found in (Bartolo et al. 2004). Testing for non-Gaussianity has then become the basic 
tool to discriminate among different models for the production of energy-density perturbations. It has become common 
practice to quantify the amount of non-Gaussianity by means of the dimensionless non-linearity parameter /nl (see, e.g. 
| |Komatsu and Spergel 2001| |), setting the strength of quadratic non-linearities in an expansion of the large-scale gravitational 
potential 4? (conventionally defined so that the temperature anisotropy is AT/T = — in the Sachs-Wolfe limit) in terms 
of a Gaussian random field (Hg, namely 

$(x) = $g(x) +/nl'4>g(x) (1) 

(up to a constant offset, which only affects the monopole contribution). Detailed calculations of the non-linearity parameter 
/nl during and after inflation (see | |Bartolo, Matarrese and Riotto 2004|IBartolo et al. 2004(1 1 have shown that it unavoidably 
contains an angle-dependent part, whose role could be extremely important to look for specific signatures of inflationary 
non-Gaussianity as recently shown in ( [Liguori et al. 2005| |. In what follows, however, we will follow the common practice of 
taking /nl as a constant parameter. The upper limits on the estimated value of /nl have become more and more stringent 
as the sensitivity of CMB experiments has improved. With MAXIMA data, IlSa.ntos et a,1. 20(^1 put a limit of |/nl| < 950, at 
1 a level. Using GOBE data, IIKomatsn et al. 200211 (with the bispectrum) and I jGayon et al. 200^ (with Spherical Mexican 
Hat Wavelets (SMHW)) found |/nl| < 1500 and |/nl| < 1100 respectively, at 1 cr level; these intervals were shrinked to —58 < 
/nl < 134 at 2(7 level in the first release of WMAP see lIKomatsu et al. 200311 . Using SMHW, i jMukherjee and Wang 2004[ 
found /nl = 50 ± 80 at 1 sigma and /nl < 220 at 2(7 level; in IIGabella et al. 200511 we constrained /nl = —5 ± 175 at 
2(7 level, combining the local curvature and spherical wavelets. Very recently the constraints on /nl have been improved by 
lIGreminelli et al. 200511 . who find — 27 < /nl < 121 at 2a level. In i jGaztanaga and Wagg 2003| l very stringent limits are found 
but a direct comparison with the other methods is unfeasible because of discrepancies in the non-Gaussian models adopted. 
In this paper we focus on functionals of the normalized bispectrum; the latter has been considered by many authors in the 
literature, including ( [Komatsu and Spergel 2001[ ; more recently, (IBabich 2005t has discussed conditions under which the 
bispectrum is the optimal estimator of primordial non-Gaussianity and ( [Babich and Zaldarriaga 2004[ showed that tighter 
constraints are expected from a joint analysis of temperature and polarization data. 

Here, we implement some procedures which were proposed in the statistical literature (IMarinncci 2005(1 . In that paper, a full 
analytic derivation is provided for the test behaviour in the presence of an ideal experiment; the behaviour in the presence 
of non-Gaussianity is also discussed. Here, we investigate the properties of these tests under a realistic experimental setting, 
using both simulations and WMAP data. The plan of this paper is as follows: in Section|^we describe the proposed procedure; 
Section El and E] describe the implementation, simulations and datasets used; Section El discusses the results of the procedure 
applied to the simulations; in SectionElwe draw some conclusions and discuss directions for future research. 


2 THE INTEGRATED BISPECTRUM 

As well-known, the angular bispectrum defined by 

As shown by IIHu 200 It for a statistically isotropic field it is convenient to focus on the angle-averaged bispectrum, defined by 
h h ^3 
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The distribution of the previous statistic depends on the angular power spectrum of the CMB. It is a standard practice 
to make the angular bispectrum model independent (under Gaussianity) by focussing on the normalized bispectrum, which 
we define by 


T _ ( 1 '>(b+^2 + t3)/2 

'I'sC 1 I n n n 


(4) 


The factor (—l)di+C+* 3)/2 jg usually not included in the definition of the normalized bispectrum; it corresponds, however, to 
the sign of the Wigner’s coefficients for mi = m 2 = m 3 = 0, and thus it seems natural to include it to ensure that and 

foiiiais share the same parity (see 0 ). An alternative estimator of I 11 I 2 I 3 , which uses the estimated rather than the theoretical 
bispectrum, is provided by 
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is the power spectrum of the given realization. In IMarinucci lOO^ has shown that, under Gaussian hypothesis, the two 
normalizations are equivalent. A crucial issue relates to how one can combine the information from the different multipoles 
into a single statistic. For statistically isotropic fields the bispectrum can be non-zero only for configurations where I 1 +I 2 + h 
is even and the triangle conditions hold, \li — lj\ < h < h + Ij, = li2,3. It is not difficult to see that if we avoid 

repetitions there are asymptotically L®/24 such configurations, where L denotes the highest observable multipole; it is therefore 
computationally very hard to consider the full set of bispectrum ordinates for high resolution experiments such as WMAP 
or Planck. Various solutions have been considered, see for instance JKomatsu, Spergel and Wandelt 2003| |. In this paper, our 
idea is to restrict the analysis to a subset of the bispectrum ordinates where the bulk of information on non-Gaussianity is 
condensed. 

The procedure we shall consider has been advocated in the statistical literature by IIMa.riniicci 200.5(1 : More precisely, for 
finite integers Iq > 2, K > 0 we shall consider the processes: 

lLr]-lo-K r 

J,, 3 ,K(r) = — ^ 

l=lo+K+l f '' 

where [.] denotes the integer part of a real number; 0 < r < 1 and lo is an (arbitrary but fixed) value which can be taken 
equal to 2 or 3, according to whether we wish to keep the quadrupole or not in the data. As usual, the sums are taken to be 
equal to zero when the index set is empty. AT is a fixed pooling parameter: for K = 0 we obtain the special case 
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( 6 ) 


The normalizing factors are chosen to ensure an asymptotic unit variance for all summands. In words, the strategy for 
JL;io,K{r) is to look at collapsed configurations; more precisely, for a fixed Iq we aim at maximizing the distance among 
multipoles, albeit preserving the triangle conditions h < Ij +lk- For an ideal experiment, it is shown in arinucci 2nn5t that, 
as L —> 00 , for any fixed integers Iq > 0, K > 0 


JL-,io,Kir) ^ W{r),0 < r < I, (7) 

where => denotes weak convergence and W (r) standard Brownian motion, that is, the Gaussian process with zero mean, inde¬ 
pendent increments and variance < W(r)^ >= r. The concept of weak convergence ensures that the asymptotic distribution 
can be immediately derived for any continuous functional of JL-,io,K{r)\ instance 

Pr < max JiL-,io,K{r) > x> = 2^{—x) for all x > 0., 

10<r<l ’ ’ I 

where Pr stands for the probability and '!>(.) denotes the cumulative distribution function of a standard Gaussian variate; 
these values are well-known and tabulated, and can be used to double-check the validity of Monte Garlo simulations. 

We wish now to discuss the expected power of this procedure under simplified circumstances: we shall work in the 
framework of an ideal experiment and a pure Sachs-Wolfe model. The Monte Garlo evidence presented in the next section 
suggests, however, that our conclusions have a much more general validity. From 0 we know that, approximately 


Var {JiL-,io,Kir)} ~ r,forlargeI/; 
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also, it is well-known that the Sachs-Wolfe bispectrum can be approximated by 


-BiPais = G/Nchiiiai, 
where G is a positive constant, 
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and lower order terms are neglected. We shall take Ci oc Z““ (for some positive constant a > 2). For simplicity, let us assume 
that the normalizing angular power spectrum is known a priori; without loss of generality, we take K = Q. Using expressions 
8.1.2.12 and 8.5.2.32 in llVarsha.lovich et al. 19881 . it can be shown that, for fixed Zo > 2, 

( 

for some C > 0 which depends on Zq but not on Z. Then we have easily that 

< Ji(r) >oc ^ ^ \/Z.oc/ nlB. (9) 

On the other hand, by a similar argument, it is easy to show that the choice of an equilateral configuration h = h ~ h = I 
entails an asymptotic negligible power in the presence of this kind of non-Gaussianity. 

The main conclusions we can draw from this heuristic discussion are as follows: “collapsed” configurations where the 
multipoles lie on the boundary of the triangle conditions seem to have an expected power of at least an order of magnitude 
in L larger than for configurations on the main diagonal. For a pure Sachs-Wolfe model and an ideal experiment, the signal 
to noise ratio is going to increase linearly for collapsed conhgurations, whereas no improvement is expected for equilateral 
configurations. It is then natural to conjecture that very little information is lost in our procedure with respect to a full 
analysis of all bispectrum coordinates; the latter, however, is clearly unfeasible for computational reasons, and no analytic 
results are available to guide the simulations. In fact whereas calculating all the elements of the bispectrum scales as P in 
CPU time, calculating only the collapsed configurations scales as P. With WMAP data {L ~ 500) we gain a factor 10® and 
for Planck resolution it is 10^ times faster. The next section is devoted to check the validity of the above claims in a much 
more general and realistic setting, by means of Monte Carlo simulations. 


Zq Z Z -|- Zq 
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3 THE AVERAGED INTEGRATED BISPECTRUM OF SIMULATED NON-GAUSSIAN MAPS 

We have generated a set of 200 non-Caussian simulations to test the power of the collapsed conhgurations described in the 
previous sections. Note that the above deductions were done for a model with fluctuations coming entirely from the Sachs-Wolf 
term with no radiative transfer involved. As the analytical treatment becomes much more complex with the radiative transfer 
function included, we will show that the above results are still valid using simulated non-Caussian maps. First of all, we will 
show that for the ’Collapsed Conhguration Bispectrum’ (CCB) signal to noise is increasing for increasing multipoles, thus 
making the detection probability monotonically rising with increasing angular resolution. Second, we will show that the power 
of the CCB is falling with increasing Zq and that all the power can be extracted using only the first few Zq. Finally, we will 
show that the often used diagonal conhguration of the bispectrum has minimal power compared to CCB. 

The non-Caussian simulated maps were generated using the method of | |Liguori, Matarrese and Moscardini 2003| l. We pro¬ 
duced a set of 200 maps (unfortunately, this is very CPU demanding limiting the number of maps which can be produced), the 
highest multipole being L = 2000 with a power spectrum similar to the best ht WMAP power spectrum IIITiushaw et al. 200311 . 
The non-Caussian part of the bispectrum can be obtain by making the ensemble average over 200 simulations of 

Lo i 

^ 0=2 + l 

where L is the maximum multipole used dependent on the resolution of the experiment and Lq is the maximum value 
of Iq included. Here NG means that we have only considered combinations of the type in calculating 

the bispectrum. Since the non-Caussian are many orders of magnitude smaller than the corresponding gaussian 
combinations with more factors will be negligible, also confirmed by our simulations. Note that in order to normalize the 
Hem, we have assumed that we are able to estimate the power spectrum well and we use the correct ensemble averaged power 
spectrum. In the case of cut sky and noise, this is taken into account in the normalizing power spectrum by Ce = CeBf / fsky+Nt 
where Be is the beam, Ne is the noise power spectrum and /sky is the sky fraction. 

In figure^ we show the Ji,io{L = 2000) for different Iq. As expected, the non-Caussian term is monotonically increasing 
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Figure 1. The non-Gaussian part of the integrated bispectrum (see eq.|2J averaged over 200 Monte Carlo simulations and normalized by 
the Gaussian standard deviation. Different values for Iq (with constant step A£o = 1) are shown, the lines are descending for increasing 
Iq, the highest being £o = 2, the lowest being £o = 16. 


and the contribution from high lo is decreasing. We also estimated the diagonal bispectrum from the same simulations and 
found it to be oscillating around 0, being at least three order of magnitude smaller than CCB. 


4 SIMULATED MAPS AND THE WMAP DATA 

In the following, we will demonstrate the power of CCB on simulated maps with noise and beam specifications being those 
of the WMAP experiment IIBenuett et al. 200.11 (all data and templates are publicly available on the LAMBDA website*^). 
Finally, we will also analyse these data and compare our results to those of tKnmat.su et al. 2nn,‘lt in which all elements of 
the bispectrum are used. We will simulate the three CMB dominated WMAP channels, the Q(41 GHz), V(61 GHz) and 
W(94GHz) bands, convolving with the corresponding beams and adding Gaussian noise. We will also use the KpO galactic 
cut. Note that the galactic cut may shift the optimal configurations of the bispectrum and as we will show later, using a 
galactic cut seems to render the CGB less optimal. All analysis will be performed on the noise-weighted linear combination 
of the Q, V and W channels, co-added according to JBennett et al. 200311 

Ti = {T^wq + Tf wv + ww)/{wq + wv + ww), 

where is the temperature in pixel i for channel X, and the weights are given as wx = l/crx where ax is the average noise 
variance for channel X. We will simulate 200 Gaussian realizations of CMB with the corresponding 200 non-Gaussian maps 
| |Liguori, Matarrese and Moscardini 2003^ as well as 200 independent Gaussian simulations and noise for all 400 maps to be 
analysed. 


5 ESTIMATING /nl IN MONTE CARLO SIMULATIONS 

The scope of this section is to find the error bars on /nl for the simulated WMAP data described in the previous section. 
These error bars will be compared to the error bars obtained by IIKomatsu et al. 200311 using all elements of the bispectrum 
as a test of the optimality of the GCB. We will apply the full estimation procedure to maps with and without galactic cut. 


^ http://lambda.gsfc.nasa.gov/ 
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Figure 2. Cross sections of the normalized correlation matrix. The plot shows ^CigigCfiii for some values of (q. The correlation 

matrix was obtained by 200 Gaussian simulations of the integrated bispectrum C(£i = {JiJgi). 


checking in this way whether introducing a sky cut causes loss of optimality (thus that the optimal configurations are shifted 
away from the CCB by the non-ortonormality of the spherical harmonic functions on the cut sky). Finally we will estimate 
/nl using the real data. 

In order to estimate /nl from the CCB of the 200 simulated maps, we will minimize the defined by 
X^(/nl) = d(/NL)^C“^d(/NL), 

with respect to /nl. Here, the elements of the data vector are defined by de = Jt — /nl * ^nd the correlation matrix 

Cei' = (dedii) — {dt){di') with Ji being the ’observed’ CCB. The 200 pure non-Gaussian maps are used to obtain and 

the 200 independent Gaussian simulations are used for obtaining the correlation matrix. Strictly speaking, the correlation 
matrix also depends on /nl but for realistic values of /nl (< 100) this dependence is negligible. The normalized correlation 
matrix defined as Cu' = Cw ! VCelCFii is shown in figure The correlations between neighbouring multipoles is so strong 
that the matrix is numerically ill-defined. A binning procedure is necessary to enable the matrix to be inverted. Note from 
the figure that the long-range correlations are stronger for larger multipoles, thus a tighter binning for the lower multipoles 
will be allowed. We define the bin-size for a given multipole I following this procedure; 

• define a limit a < 1 

• start with 1 = 2 and find for which I the value of the normalized correlation matrix C 2 t has fallen below the limit a. This 
defines the next bin. 

• Starting with the obtained multipole i' of the next bin, find for which multipole C^i has fallen below the limit a. 

• repeat the above procedure until the highest multipole L has been reached. Check if the correlation matrix with this 
binning gives a numerically well defined correlation matrix. If not, repeat the above procedure with a lower limit a otherwise, 
the binning procedure is finished and the final binning has been obtained. 

Following this procedure, we obtained the binning shown in the panels of figure |3| for the different cases which we will 
describe in the following. Further, we estimate /nl in the 200 Gaussian maps and obtain in this way the frequentist error 
bars. We have checked that the we obtain unbias estimates of /nl and that error bars of maps with non-zero /nl for realistic 
values (/nl < 100) do not deviate significantly from the value obtained on Gaussian realizations. 

As a first test of this procedure, we estimated /nl for an ideal experiment with no noise and no sky cut, L = 2000 and 
Lo = 16. In this case we obtained A/nl = 8 at Ict being consistent with ( [Komatsu and Spergel 2001| who found A/nl = 3 
for an ideal experiment with L = 3000 using all elements of the bispectrum. 
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Figure 3. Left panel: Multipole binning for the analysis of the high resolution simulations. On the x-axis we show the bin number 
and on the y-axis the corresponding last multipole of the bin. The crosses show each 16th bin to give an impression of the bin density 
per multipole. The limit a on the correlation matrix used to construct this binning (see the text) was a = 0.998. Right panel: the same, 
but for the WMAP simulations. The crosses are now shown for all bins. Lower line is for the case with WMAP noise but no galactic cut 
(a = 0.985, see the text), upper line is for the case with both noise and a KpO cut (a = 0.946) 


Second, we applied the same procedure on the maps with WMAP beam and noise added (but no galactic cut) and find 
^/nl = 40 which is to be compared with A/nl = 48 obtained by (IKomatsu et al. 2003t for the WMAP data. In the latter, 
a galactic cut of 25% (KpO) was applied. 

Finally, for the case including the KpO sky cut, we obtain A/nl = 80. Note that in this case the error has increased 
drastically and much more than one could expect from a simple 25% loss of data. This can be understood, looking at the 
right panel of figure showing the binning with and without the mask. The bins are much denser without the mask, showing 
that the mask is introducing high correlations between multipoles. Thus the information is spread to other multipoles and 
the collapsed configuration is no longer optimal. This suggests to include refilling procedure in order to restore orthogonality 
(i.e. one could use the Gibbs sampling approach IIEriksen et al. 200411 ) and just be limited by the sampling variance, this will 
be explored in a future work. 

We have applied the procedure to the WMAP data and obtain an estimate /nl = 0 when sampling /nl space in a grid 
of 10. In figure 0 we show the Gaussian standard deviation of Je together with Ji for the WMAP data. We also show the 
theoretical Ji for a set of /nl values. On the right panel of figure |11 we show the yf around its minimum, showing that the 
Bayesian error bars indicate error bars of A/nl about 80 at la and 160 at 2a in agreement with the frequentist analysis. 


6 CONCLUSIONS 

The bispectrum is one of the most common statistics to test for non-Gaussianity in CMB data. In particular, for estimating 
the non-linear coupling constant in primordial non-Gaussian models, the bispectrum has proven to produce very stringent 
limits. The drawback of the bispectrum is that the number of elements scale with the maximum multipole L as and 
computational time scales as L®. 

In this paper we have introduced a new procedure to test for non-Gaussianity on GMB data; our approach is based 
on an integrated form of the bispectrum, where a phase factor is introduced and the focus is narrowed on nearly-collapsed 
configurations. Our approach is computationally very convenient as it scales only as L® and presents the added bonus to allow 
for explicit analytic results under idealized circumstances, in both Gaussian and non-Gaussian settings. The power properties 
are shown to be encouraging by Monte Carlo experiments: indeed some simple calculations suggest that the non-Gaussian 
signal grows linearly with the experiment resolution. 

By comparing to limits on /nl published in the literature based on all elements of the bispectrum, we have shown that the 
collapsed configuration bispectrum does indeed seem to produce comparable near optimal results. For the ideal experiment 
we can constrain /nl < 8 using L = 2000 and for a WMAP like experiment we find /nl < 40 when a galactic cut is not 
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Figure 4. Left panel: The integrated bispectrum (see eq. |3 averaged over 200 Monte Carlo simulations with WMAP noise and KpO 
mask. The solid line represents the mean of 200 Gaussian simulations, dotted (dashed) lines represent non-Gaussian simulations for 
different negative (positive) values of /nl from -500 to 500. The shaded area indicates the 1 a confidence level taken from the Gaussian 
simulations. The diamonds show the result of WMAP data. Right panel: The of WMAP data as a function of /nl- /nl is estimated 
to be 0 it 80 and 0 ± 160 at Icr and 2t7 level respectively. 


introduced. However, it turns out that a galactic cut does destroy the nice properties of the collapsed configuration. The 
multipoles get strongly coupled and the constrains on /nl becomes bigger that what one would expect from a pure increase 
in sampling variance. For this reason, when the galactic cut is introduced, we obtain —80 < /nl < 80 at Icr applied to the 
WMAP data. Clearly, data from future CMB experiments like Planck will still need to be analysed with a galactic cut as it 
will be impossible to completely eliminate the galactic plane. For that reason, a way of dealing with the non-optimality of 
the collapsed configuration bispectrum will be necessary. One possible solution to this would be by some refilling procedure 
(or one could use the already existing Gibbs sampling technique IIFrikseu et al. 20041 1 which could restore orthogonality of 
the spherical harmonics. In this work, we did not make any prediction of the constrains on /nl which can be achieved by 
the Planck experiment. As we are still unable to deal optimally with cut sky, this would not produce a precise limit and is 
therefore postponed to future work. 

Finally, we note that for an ideal experiment it seems feasible to achieve the bounds predicted in ( jKomatsu and Spergel 2001} , 
despite the fact that we are using here only L, rather than H, bispectrum configurations. On the other hand, the presence of 
gaps greatly deteriorated the performance of our procedure. As the collapsed bispectrum seems potentially a very promising 
technique for high resolution data, we view the gap handling issue as a research priority in our future work. 
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